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1.  Introduotlon 


We  will  show  how  to  generate  normal  remdom  variables  very 
rapidly  In  a  oomputer  -  for  example,  at  the  rate  of  10,000  -  15,000 
per  second  In  the  IBM  7090.  The  method  is  suitable  for  any  computer. 
It  has  evolved  from  a  number  of  procedures  we  have  considered  In 
the  past,  [1]  -  [3].  The  Incorporation  of  successive  Improvements 
has  led  to  a  procedure  which  1)  Is  fairly  esisy  to  program,  2)  requires 
little  storage,  300  -  400  constants,  3)  is  very  fast  -  It  takes 
about  as  long  to  generate  the  normal  x  aa  the  \inlform  u  from 
%dilch  It  comes,  and  4)  le  completely  eiccurate,  in  the  sense  that  In 
theory  the  procedure  returns  a  random  variable  with  exactly  the 
required  distribution;  In  practice  the  result  Is  an  approximation 
Influenced  only  by  the  capacity  (word  length)  of  the  computer. 

In  short,  our  method  Is  much  faster  than  any  we  have  heard  of, 
and  Is  completely  accurate.  We  recommend  that  It  be  used  as  the 
basic  normal  random  variable  generator  in  any  computer  Installation. 
Comparisons  of  methods,  times,  storage  requirements,  etc.,  are 
given  In  Section  4. 

2.  The  procedure  for  decimal  computers 

We  will  go  Into  detail  for  the  decimal  ceuse.  The  method  for  the 
binary  case  Is  similar  and  the  necessary  modifications  will  be  given 


2 


in  Seotlon  3.  Suppose  we  have  a  procedure  for  generating  vmlform  [0»1] 
random  variables  u  with,  say,  8  decimal  places.  Then  probably 
the  fastest  way  to  generate  a  variable  x  with  absolutely  continuous 

g 

distribution  F  is  to  assign  a  value  of  x  for  each  of  the  10 

possible  u's  according  to  the  relation  x  =  F"'^(u)  .  The  time 

for  generating  x  would  then  be  the  time  to  generate  u,  look  up 

the  value  stored  in  the  location  associated  with  that  u,  and  bring 

it  to  the  required  place.  Unfortunately,  this  fastest  procedure 
o 

would  require  10  storage  locations.  Our  procedure  approaches 
this  ultimate  in  speed  in  the  following  way:  96  -  97^  of  the  time, 
we  use  the  first  few  digits  of  u  to  locate  one  of  only  a  few 
hundred  stored  values,  then  we  add  the  last  digits  of  u  to  that 
stored  value.  The  other  3  -  .U%  of  the  time,  we  do  what  is  necessary 
to  make  the  resulting  mixture  come  out  right.  The  average  time  for 
the  entire  procedure  still  remains  close  to  the  ultimate  which  can 

g 

be  attained  by  using  10  storage  locations.  The  compact  storage 
procedure  which  enables  us  to  conseirve  storage  space  is  described  in 

[4]. 


Now  for  details  of  the  method.  Assume  we  have  a  procedure  for 
producing  Independent  xiniform  [0,1]  random  variables 

u,U2,U2,u^, .... 

The  problem  is  to  generate  a  normal  random  variable  X  in  terms  of 
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i 

1 


Figure  1 
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the  u'e.  Let  g  be  the  absolute  normal  density t 

_^2 

(1)  g(*)  =  ^—r  *  i  0 

*y2TT 

The  procedures  below  provide  a  random  variable  with  density  g)  a 
random  +  must  be  attached  at  some  convenient  point  in  the  progrcun, 
or  else,  by  doubling  the  storage  requirements,  the  procedure  may  be 
modified  to  handle  the  standard  normal  density  on  - <»  <  x  <  ». 

We  write  g  as  a  mixture  of  three  densities,  as  pictured  in 
Figure  It 

(2)  g(x)  =  .9578g^(x)  +  .0395g2(x)  +  .0O27g3(x) . 

The  probabilities  displayed  in  (2)  have  been  rounded  to  4  places | 
they  are  actually  used  to  10  places  in  the  program.  With  probability 
.9578. ••  we  generate  a  random  variable  with  density  g^,  with 
probability  .0393. ••  we  generate  a  random  variable  with  density 
g^,  and  with  probability  .0027...  we  generate  a  random  variable 
with  density  g^.  The  time  for  g^  is  very  short,  for  g2  short, 
and  for  g^,  long. 

We  may  produce  a  variable  with  density  g^  in  the  form  z  +  .lu, 
where  z  is  discrete  and  u  uniform  [0,1].  The  oaapaot  storage 
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method  described  in  [4]  provides  z  qulokl7,  with  modest  storage 
requirement.  The  altitudes  of  the  rectangles  in  are  truncated 
to  three  decimal  places  for  use  In  that  method,  the  remaining 
portion  Is  lumped  with  the  toothlike  region  abo/e  the  rectangle 
and  dealt  with  directly  on  those  Infrequent  occasions  when  It  Is 
required. 

A  variable  with  density  g^,  required  about  4  percent  of 
the  time,  Is  produced  as  follows:  One  of  the  "teeth"  from  g^ 

Is  selected  with  appropriate  probability,  then  a  random  variable 
with  the  nearly  linear  density  given  by  that  tooth  Is  generated. 

The  procedure,  described  In  [3],  Is  summarized  as  follows: 

To  generate  a  random  variable  X  with  a  nearly  linear 
density  function  g(x) ,  s  ^  x  ^  s  o 

such  as  this  or  this 


enclose  g(x)  within  two  parallel  lines. 
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Ilk*  this 


or  this. 


Them 

1.  Choose  Independent  \inlform  [0,l]  rtmdom  variable*  u  and  v* 
2*  If  niax(u,v)  <  a/b  put  X  =  s  +  c  inln(u,v) . 

3«  If  not,  test!  b|u  -  v|  ^  g[8  +  c  niin(u,v)]?  if  yes,  put 
X  =  a  +  c  oln(u,v)}  If  no,  go  to  step  1  and  try  again. 

In  our  application  of  this  method,  we  have  c  >  .1  and  the 
teet  in  step  3  may  be  put  in  the  form 

k|u  -  v|  <  -  1,  where  w  =  -  |{[b  +  .1  mln(u,v)]^  -  [s  +  .1]^} 

and  k  =  .1(8  +  .1)  for  a  =  0, .1, .2, . . . , .9  and 
k  =  _  2  for  a  =  1.0,1. 1,1. 2,..., 2.9. 

The  ratios  a/b  are  stored  in  memory  locations  230-259.  Moat 
of  them  are  close  to  1,  so  that  step  2  of  the  nearly  linear  technique 
will  provide  X  most  of  the  tlmei  only  occaaloncdly  will  an  exponen- 
ticd  subroutine  be  required. 
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Density  comes  from  the  tall  of  g.  It  will  provide  our 
random  variable  about  1  time  in  400.  We  then  need  an  absolute  normal 
variable,  |x|,  conditioned  by  |x|  >  3.  We  get  it  by  choosing  one 
of  a  pair  of  absolute  normal  variables  lying  in  the  first 

quadrant  but  outside  the  square  t 


We  generate  by  choosing  a  normal  point  outside  the  quarter 

circle,  rejecting  it  if  it  lies  inside  the  square.  Thus,  if 

***•  l-“*®P«“dent,  uniform  [0,1],  conditioned  by  uj  +  Ug  ^  1, 

ve  put 


=  u  r-2±2sLii 
“ll  2  2^ 


“l*“2 


'■  2  2 
Ui^Ug 


where  w  has  the  exponential  distribution.  We  may  provide  w  in  the 
2  2  2  2 

form  -  ln(u--fu2) ,  since  u^+u^  is  luiiform  [0,1]  and  is  independent 
^  (see  [5]).  Then  if  (x^,X2)  lies  outside  the  square,  we  take 
whichever  of  the  coordinates  that  is  ^3.  If  both  are  ^3,  we  might 
possibly  store  one  for  future  use,  but  the  relative  frequency  for  both 
^  3,  about  1  in  800,  is  not  enough  to  justify  asking  for  a  possible 
stored  value.  The  measures  of  the  regions  are  given  in  the  figure,  so 
that  54/( 57+54) »  or  about  U9%  of  the  time,  a  pair  (xj^,X2)  lying 
outside  the  quarter-circle  will  lie  outside  the  square. 
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Outline  of  the  program i  The  following  outline  deecrlbes  the  procedure 
for  a  decimal  machine,  and  the  flow  chart  on  the  opposite  page  describes 
the  same  procedure  with  a  little  more  detail.  Assume  that  there  Is  a 
subroutine  for  producing  Independent  uniform  [0,1]  random  variables 
u,u^,U2,...  on  demand.  Let  the  constants  C(n)  be  stored  In  memory 
locations  n  *  0,1,2, . .  .,289.  Let  t  ^  10  .  Let  u  *  .dQ^d^d^d^d^. .  • 

be  the  first  u,  the  d's  Its  decimal  digits.  Then: 


1. 


2. 


3. 

4* 


If  00  ^ 

If  790  ^  d^djd^  <  940, 

If  .94  <  u  <  .9973002039, 

If  .9973002039  <  u  <  1, 

3^-2  ln(u?  +  u?)  1 

=  “it — r;-  2'  i* 

"l  ^"2 


put  X  =  CCdj^dj)  +  .Od^djd^... 
put  X  =  cCdj^dgd^  -  771)  +  .Od^d^d^... 
put  J  =  170  and  go  to  step  5* 
form  pairs 


*2  =  ^"2^ 


3=  -  2  l»(u=  * 


2^2 

Ui  +U2 


until  one  of  the  pair  ^>^2  '^2 

2  2 

are  conditioned  by  u^^  +  Ug  <  1. 

5.  Test*  Is  u  <  C(J)7  If  yes,  go  to  step  6.  If  no,  put  J  =  J  +  1 
and  repeat  step  5. 

6.  Testi  Is  u  <  C(J  +  30)?  If  yes,  go  to  step  8.  If  no,  go  to  step  7» 

7.  Generate  new  u  and  put  x  =  C(J  -  30)  +  6u. 

8.  Generate  new  Test:  Is  max(uj^,U2)  <  C(J  +  60)?  If  yes, 

put  X  =  C(J  -  30)  +  6  mln(uj^,U2).  If  no,  let 

w  =  -  .5[C(J  -  30)  +5  mln(Uj^,U2)]^  -  [C(J  -  30)  +  5]^ 
and  test:  Is  \vi^  -  U2I  <  C(J  +  90) (e''  -  1)?  If  yes,  put 

X  =  C(J  -  30)  +  6  mln(u^,U2).  If  no,  repeat  step  8. 


To  generate  a  noraa! 


*  .didfdi... 
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L  variable  X  In  a  decimal  o(»iiutert 


bdl 


be 
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ConstAnts  for  Loading  In  Memory  of  a  Decimal  Machine 


Location  Contents 


1 

0 

0.^ 

50 

0.5 

100 

1.6 

150 

0.5 

200 

.9422781966 

250 

.973 

1 

0.2 

51 

0.6 

101 

1.6 

151 

1.2 

201 

.9455720777 

251 

.975 

2 

0.3 

52 

0.7 

102 

1.6 

152 

1.6 

202 

.9485514463 

252 

.974 

3 

0.3 

53 

0.7 

103 

1.6 

153 

1.7 

203 

.9511653133 

253 

.978 

4 

0.3 

54 

0.7 

104 

1.7 

154 

0.3 

204 

.9549863293 

254 

.755 

5 

0.3 

55 

0.7 

105 

1.7 

155 

1.5 

205 

.9566914271 

255 

.970 

6 

0.3 

56 

0.7 

106 

1.7 

156 

2.0 

206 

.9604850173 

256 

.501 

7 

0.5 

57 

0.8 

107 

1.8 

157 

1.8 

207 

.96380a343 

257 

.971 

8 

0.6 

58 

0.8 

108 

1.9 

158 

2.2 

208 

.9665717757 

258 

.968 

9 

0.6 

59 

0.9 

109 

1.9 

159 

0.2 

209 

.9689169701 

259 

.967 

10 

0.6 

60 

0.9 

no 

1.9 

160 

2.5 

210 

.9712916782 

260 

12.5 

11 

0.6 

61 

0.9 

111 

1.9 

161 

2.3 

211 

.9742012516 

261 

8.20523339 

12 

0.6 

62 

0.9 

112 

1.9 

162 

2.4 

212 

.9761328124 

262 

6.91865398 

13 

0.8 

63 

1.0 

113 

1.9 

163 

2.1 

213 

.9784228835 

263 

20.0 

14 

0.8 

64 

1.0 

114 

1.9 

164 

0.1 

214 

.9805795259 

264 

,9.03255791 

15 

0.8 

65 

1.1 

115 

1.9 

165 

2.7 

215 

.9830652062 

265 

4.64444483 

16 

1.0 

66 

1.1 

116 

2.0 

166 

0.0 

216 

.9842240767 

266 

6.40863082 

17 

1.0 

67 

i.i 

117 

2.0 

167 

2.6 

217 

.9863251515 

267 

10.0 

18 

1.5 

68 

1.2 

118 

2.0 

168 

2.8 

218 

.987ia5820 

268 

11.11111111 

19 

0.0 

69 

1.2 

119 

2.0 

169 

2.9 

219 

.9888328519 

269 

14.2857142 

20 

0.0 

70 

1.2 

120 

2.0 

170 

.9432165072 

220 

.9894907759 

270 

16.666^666 

21 

0.0 

71 

1.3 

121 

2.0 

171 

.9464092887 

221 

.9907816118 

271 

7.510a395 

22 

0.0 

72 

1.3 

122 

2.0 

172 

.9494969394 

222 

.9917305989 

272 

5.57434982 

23 

0.0 

73 

1.4 

123 

2.1 

173 

.9525783787 

223 

.9930632865 

273 

5.22886160 

24 

0.0 

74 

1.4 

124 

2.1 

174 

.9555567647 

224 

.993813a06 

274 

25.0 

25 

0.0 

75 

1.5 

125 

2.1 

175 

.9584896204 

225 

.9942625460 

275 

5.96452440 

26 

0.1 

76 

1.6 

126 

2.1 

176 

.9613885364 

226 

.9951108014 

276 

4.39512015 

27 

0.1 

77 

1.7 

127 

2.1 

177 

.9641982792 

227 

.9958055523 

277 

4.92081328 

28 

0.1 

78 

1.8 

128 

2.1 

178 

.9667888257 

228 

.9960778669 

278 

3.96317864 

29 

0.1 

79 

0.0 

129 

2.2 

179 

.9693677568 

229 

.996a38342 

279 

33.33333333 

30 

0.1 

80 

0.4 

130 

2.2 

180 

.9719365988 

230 

.973 

280 

3.44279563 

31 

0.1 

81 

0.4 

131 

2.2 

181 

.97U749700 

231 

.996 

281 

3.77488448 

32 

0.1 

82 

0.7 

132 

2.2 

182 

.9769426279 

232 

.992 

282 

3.60202892 

33 

0.2 

83 

0.9 

133 

2.3 

183 

.9792129152 

233 

.920 

283 

4.16906566 

34 

0.2 

84 

0.9 

134 

2.3 

184 

.9812335540 

234 

.998 

284 

50.0 

35 

0.2 

85 

0.9 

135 

2.3 

185 

.9832493731 

235 

.982 

285 

3.15925147 

36 

0.2 

86 

1.1 

136 

2.4 

186 

.9850207959 

236 

.990 

286 

100.0 

37 

0.2 

87 

1.1 

137 

2.4 

187 

.9864483145 

237 

.996 

287 

3.29564243 

38 

0.3 

88 

1.1 

138 

2.5 

188 

.9878069895 

238 

.985 

288 

3.03248984 

39 

0.3 

89 

1.1 

139 

2.6 

189 

.98911oa50 

239 

.959 

289 

2.91437825 

40 

0.4 

90 

1.3 

140 

0.7 

190 

.9902073697 

240 

.942 

a 

0.4 

91 

1.3 

141 

1.1 

191 

.9912605179 

2a 

.994 

42 

0.4 

92 

1.3 

142 

1.3 

192 

.9922362590 

242 

.986 

43 

0.4 

93 

1.3 

143 

0.4 

193 

.9931582051 

243 

.985 

U 

0.4 

94 

1.3 

144 

1.0 

194 

.9940219494 

244 

.890 

45 

0.4 

95 

1.3 

145 

1.9 

195 

.9948456363 

245 

.988 

46 

0.5 

96 

1.4 

146 

1.4 

196 

.9955013109 

246 

.980 

47 

0.5 

97 

1.4 

147 

0.9 

197 

.9958897393 

247 

.983 

48 

0.5 

98 

1.6 

148 

0.8 

198 

.9962683734 

248 

.977 

49 

0.5 

99 

1.6 

149 

0.6 

199 

.9973002039 

249 

.843 

i 
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3»  The  procedure  for  binary  computers 

The  methcxl  for  binary  machines  Is  much  the  same  as  for  the  decimal 
case}  the  constants  change,  of  course,  and  the  first  four  octal  digits 
of  u  are  used  to  locate  a  discrete  variable,  rather  than  the  first 
three  in  the  decimal  case.  The  binary  program  requires  456  storage 
locations  (456  base  10  =  707  base  8) .  Here  Is  the  outline  for  a  blpary 
machine)  all  numbers  In  the  outline  are  octal  except  the  number  of  the  steps 
Let  u  =  .b^bgb^b^b^.. .  be  a  uniform  [0,1]  random  variable,  the 
b's  Its  octal  digits.  Let  the  constants  C(n)  be  stored  in. memory 


locations  n  =  0,1,2, ... ,707.  Lot 

h  = 

2"^. 

1.  If 

00  <  b^bj  <  54, 

pwt 

X  =  C(b^b2)  +  (.bjb^...)6 

2.  If 

540  ^  bj^b2b3  <  733, 

put 

X  =  C(b^b2b2  -  506)  +  (.bjb^...)6 

3.  If 

7330  i  <  7571, 

put 

X  =  C(b^b2b2b^  -  7161)  +  (.bjb^...)fc 

4.  If 

.7571  ^  u  <  .776474207403, 

put 

J  =  330  and  go  to  step  6. 

5.  If 

.776474207403  ^  n  <  1, 

form 

pairs  Xj^,X2S 

*1 

3^  -  2  In  u-  1 

"l  ''2 

3^  -  2  In  u,  1 

"l*”2 

lutll  one  of  the  pair  ^)3C2  Is 

^3, 

and  let  that  be  x)  and  U2 

2  2 

are  conditioned  by  u^  +  U2<1. 

6.  Test!  Is  u  <  C(J)?  If  yes,  go  to  step  7.  If  no,  put  J  =  J  +  1 
and  repeat  step  6, 

7.  Test!  Is  u  <  C(J  +  60)?  If  yes,  go  to  step  9.  If  no,  go  to  step  8. 

8.  Generate  nev  u  and  put  x  =  C(j  -  60)  +  &u. 

9.  Generate  new  '*'*®^*  <  C(J  +  140)?  If  yes, 

put  X  =  C(J  -  60)  +  t  min(uj^,U2).  If  no,  let 

w  =  -  .5[C(J  -  60)  +  6  min(uj^,U2)]^  -  [C(J  -  60)  +  6]^ 

and  tests  Is  |uj^  -  Ujl  <  C(j  +  220) (e''  -  1)?  If  yes,  put 

X  s  C(J  -  60)  +  6  inln(u^,U2)«  If  no,  repeat  step  9* 
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To  generate  a  normal  variable  X  In  a  binary  computer! 
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4>  Rsmarka 

The  fast  parts  of  the  program,  steps  1-2  for  the  decimal  and 
1-3  for  the  binary,  are  used  most  of  the  time.  They  should  be  vxltten 
in  machine  language,  with  care.  A  procedure  for  assigning  a  random 
+  must  be  incorporated,  too.  The  remaining  steps  can  be  written  in 
FORTRAN  or  some  such  language  with  little  effect  on  the  average  running 
time. 


We  have  written  programs  based  on  this  procedure  for  the  IBM  1620 
(decimal) ,  and  the  IBM  7090  (binary)  machines .  The  programs  are 
written  as  standard  FORTRAN  function  subprograms  with  the  necessary 
linkages,  setting  index  registers,  returning  X  in  normalised  floating 
form,  etc.  In  the  IBM  7090,  each  call  to  the  subroutine  for  a  normal 
variable  requires  about  50  cycles,  which  gives  a  rate  of  about  10,000 
per  second.  If  the  procedure  is  Incorporated  in  a  larger  program  and 
not  written  as  a  standard  subroutine,  then  rates  of  about  15,000  per 
second  are  possible.  The  standard  subroutine  for  the  IBM  7090  requires 
about  1,300  storage  locations,  including  the  space  for  the  constants 
and  for  the  instructions. 

The  IBM  1620  decimal  machine  is  quite  a  bit  slower.  The  above 
procedure,  written  as  a  standard  FORTRAN  subprogram  for  that  machine, 
with  necessary  linkages,  returning  X  in  normalized  floating  point 
form,  etc.,  takes  about  20  milliseconds.  The  constants  and  the 
instructions  require  about  4,000  core  storage  positions. 
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